week 7: multilevel models
multilevel adventures
Code
set.seed(9)
sample_id = sample(unique(d$id), replace=F, size = 20)
distinct(d, id, female, day) %>%
filter(id %in% sample_id) %>%
add_epred_draws(m3) %>%
mean_qi() %>%
ggplot( aes(x=day, y=.epred, group = as.factor(id), color=as.factor(female))) +
geom_line() +
scale_color_manual(values=c("#1c5253" , "#e07a5f")) +
labs(x="day",y="PA.std", color = "female") +
facet_wrap(~female) +
theme(legend.position = "top")
more than one type of cluster
McElreath doesn’t cover this in his video lecture, but this is from the textbook and worth discussing.
Data from Silk et al. (2005)
Zooming in on just the actor and block effects. (Remember, these are differences from the weighted average.)
Code
as_draws_df(m4) %>%
select(starts_with("sd")) %>%
pivot_longer(everything()) %>%
ggplot(aes(x = value, fill = name)) +
geom_density(linewidth = 0, alpha = 3/4, adjust = 2/3, show.legend = F) +
annotate(geom = "text", x = 0.67, y = 2, label = "block", color = "#5e8485") +
annotate(geom = "text", x = 2.725, y = 0.5, label = "actor", color = "#0f393a") +
scale_fill_manual(values = c("#0f393a", "#5e8485")) +
scale_y_continuous(NULL, breaks = NULL) +
ggtitle(expression(sigma["group"])) +
coord_cartesian(xlim = c(0, 4))
Let’s look at the predictions by actor and block to confirm.
Code
d2 %>% distinct(actor, block, treatment) %>%
add_epred_draws(m4) %>%
mutate(treatment=factor(treatment,
levels=as.character(1:4),
labels=t_labels)) %>%
group_by(actor, treatment) %>%
mean_qi(.epred) %>%
ggplot( aes(x=treatment, y=.epred, group=1) ) +
geom_point() +
geom_line() +
labs(x=NULL, y="p(pull left)", title="by actor") +
facet_wrap(~actor)
Code
d2 %>% distinct(actor, block, treatment) %>%
add_epred_draws(m4) %>%
mutate(treatment=factor(treatment,
levels=as.character(1:4),
labels=t_labels)) %>%
group_by(block, treatment) %>%
mean_qi(.epred) %>%
ggplot( aes(x=treatment, y=.epred, group=1) ) +
geom_point() +
geom_line() +
labs(x=NULL, y="p(pull left)", title="by block") +
facet_wrap(~block)
LKJ correlation prior
Code
# examples
rlkj_1 = rethinking::rlkjcorr(1e4, K=2, eta=1)
rlkj_2 = rethinking::rlkjcorr(1e4, K=2, eta=2)
rlkj_4 = rethinking::rlkjcorr(1e4, K=2, eta=4)
rlkj_6 = rethinking::rlkjcorr(1e4, K=2, eta=6)
data.frame(rlkj_1= rlkj_1[,1,2],
rlkj_2= rlkj_2[,1,2],
rlkj_4= rlkj_4[,1,2],
rlkj_6= rlkj_6[,1,2]) %>%
ggplot() +
geom_density(aes(x=rlkj_1, color = "1"), alpha=.3) +
geom_density(aes(x=rlkj_2, color = "2"), alpha=.3) +
geom_density(aes(x=rlkj_4, color = "4"), alpha=.3) +
geom_density(aes(x=rlkj_6, color = "6"), alpha=.3) +
labs(x="R", color="eta") +
theme(legend.position = "top")
More about stat_ellipse here.
We can use the slopes and intercepts to calculate the expected morning and afternoon wait times for each cafe.
Code
cafe_params %>%
mutate(
morning = intercepts,
afternoon = intercepts + slopes
) %>%
ggplot( aes(x=morning, y=afternoon) ) +
mapply(function(level) {
stat_ellipse(geom = "polygon", type = "norm",
linewidth = 0, alpha = .1, fill = "#1c5253",
level = level)
},
# enter the levels here
level = c(1:9 / 10, .99)) +
geom_point(size = 2)+
labs(x="morning wait time",
y="afternoon wait time")
Or, we can use epred_draws.
Code
d3 %>%
distinct(cafe, afternoon) %>%
add_epred_draws(m5) %>%
group_by(cafe, afternoon) %>%
summarise(wait = mean(.epred), .groups = "drop") %>%
mutate(afternoon =ifelse(afternoon==1, "afternoon", "morning")) %>%
pivot_wider(names_from = afternoon, values_from=wait) %>%
ggplot( aes(x=morning, y=afternoon) ) +
mapply(function(level) {
stat_ellipse(geom = "polygon", type = "norm",
linewidth = 0, alpha = .1, fill = "#1c5253",
level = level)
},
# enter the levels here
level = c(1:9 / 10, .99)) +
geom_point(size = 2)+
labs(x="morning wait time",
y="afternoon wait time")What is the correlation of our intercepts and slopes?
# A tibble: 1 × 6
cor_cafe__Intercept__afternoon .lower .upper .width .point .interval
<dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 -0.506 -0.799 -0.0956 0.95 mean qi Code
post = as_draws_df(m5)
rlkj_2 = rethinking::rlkjcorr(nrow(post), K=2, eta=2)
data.frame(prior= rlkj_2[,1,2],
posterior = post$cor_cafe__Intercept__afternoon) %>%
ggplot() +
geom_density(aes(x=prior, color = "prior"), alpha=.3) +
geom_density(aes(x=posterior, color = "posterior"), alpha=.3) +
scale_color_manual(values = c("#1c5253" , "#e07a5f")) +
labs(x="R") +
theme(legend.position = "top")
multilevel people again
Let’s apply what we’ve learned to our affect data:
PA_lag is a within-person variable (Level 1), whereas female is a between-person variable (Level 2).
Code
# plot posterior density
p1 = post %>%
ggplot(aes(x = b_PA_lag)) +
stat_halfeye(fill="#1c5253") +
geom_vline( aes(xintercept=0), linetype = "dashed" ) +
labs( title="posterior density",
x="slope",
y=NULL) +
scale_y_continuous(breaks=NULL)
# posterior lines
p2 = d %>%
sample_n(2000) %>%
ggplot(aes( x=PA_lag, y=PA.std )) +
geom_point(alpha=.2, shape = 20) +
geom_abline( aes(intercept=b_Intercept, slope=b_PA_lag),
data=post[1:50, ],
alpha=.3,
color="#1c5253") +
labs( x="NA",
y="PA",
title="posterior lines")
p1 + p2 
Code
p1 = m1_sigma %>%
full_join(m6_sigma) %>%
ggplot( aes(y=sigma, fill=model )) +
stat_halfeye(alpha=.5) +
labs(
x=NULL,
y = "within-person sd",
title = "posterior distributions of\neach model") +
scale_x_continuous(breaks=NULL) +
theme(legend.position = "bottom")
p2 = data.frame(diff = m6_sigma$sigma - m1_sigma$sigma) %>%
ggplot( aes(x=diff)) +
stat_halfeye() +
labs(
y=NULL,
x = "m6-m1",
title = "posterior distribution of\ndifference in sigma") +
scale_y_continuous(breaks=NULL)
p2 + p1
individual differences
To what extent do people differ in the association of holdover PA?
To what extent do people differ in the association of holdover PA?
To what extent do people differ in the association of holdover PA?
Code
set.seed(1)
sample_id = sample(unique(d$id), replace=F, size = 6)
d %>%
filter(id %in% sample_id) %>%
ggplot( aes( x=PA_lag, y=PA.std ) ) +
geom_point(alpha=.5, shape=20) +
geom_abline( aes(intercept=b_Intercept, slope=b_PA_lag),
data = post[1:50, ],
alpha=.4,
color="#1c5253") +
facet_wrap(~id) +
labs(y="today",
x="yesterday")
Code
post2 = m6 %>% spread_draws(b_Intercept, b_PA_lag, r_id[id, term]) %>%
pivot_wider(names_from = term, values_from = r_id) %>%
mutate(intercept = b_Intercept + Intercept,
slope = b_PA_lag + PA_lag) %>%
filter(id %in% sample_id) %>%
with_groups(id, sample_n, 50)
d %>%
filter(id %in% sample_id) %>%
ggplot( aes( x=PA_lag, y=PA.std ) ) +
geom_point(alpha=.5, shape=20) +
geom_abline( aes(intercept=intercept, slope=slope, color=as.factor(id)),
data = post2,
alpha=.4) +
facet_wrap(~id) +
guides(color="none") +
labs(y="today",
x="yesterday")
Code
post2 = m6 %>% spread_draws(b_Intercept, b_PA_lag, r_id[id, term]) %>%
pivot_wider(names_from = term, values_from = r_id) %>%
mutate(intercept = b_Intercept + Intercept,
slope = b_PA_lag + PA_lag) %>%
filter(id %in% sample_id) %>%
with_groups(id, sample_n, 50)
d %>%
filter(id %in% sample_id) %>%
ggplot( aes( x=PA_lag, y=PA.std ) ) +
geom_point(alpha=.5, shape=20) +
geom_abline( aes(intercept=intercept, slope=slope, color=as.factor(id)),
data = post2,
alpha=.4) +
geom_abline( aes(intercept=b_Intercept, slope=b_PA_lag),
data = post[1:50, ],
alpha=.4) +
facet_wrap(~id) +
guides(color="none") +
labs(y="today",
x="yesterday")